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Simulating the evolution of soot mixing state with a 
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Abstract. The mixing state of soot particles in the atmosphere is of crucial importance 
for assessing their climatic impact, since it governs their chemical reactivity, cloud con- 
densation nuclei activity and radiative properties. To improve the mixing state repre- 
sentation in models, we present a new approach, the stochastic particle-resolved model 
PartMC-MOSAIC, which explicitly resolves the composition of individual particles in a 
given population of different types of aerosol particles. This approach accurately tracks 
the evolution of the mixing state of particles due to emission, dilution, condensation and 
coagulation. To make this direct stoc:hastic particle-based method practical, we imple- 
mented a new multiscale stochastic coagulation method. With this method we achieved 
optimal efficiency for applications when the coagulation kernel is highly non-uniform, as 
is the case for many realistic applications. PartMC-MOSAIC was applied to an ideal- 
ized urban plume case representative of a large urban axea to simulate the evolution of 
carbonaceous aerosols of different types due to coagulation and condensation. For this 
urban plume scenario we quantified the individual processes that contribute to the ag- 
ing of the aerosol distribution, illustrating the capabilities of our modeling approach. The 
results showed for the first time the multidimensional structure of particle composition, 
which is usually lost in internally-mixed sectional or modal aerosol models. 
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1. Introduction 

Soot particles are an important constituent of the atmo- 
spheric aerosol, since they participate in troposphcric chem- 
istry [Saathoff et al, 2001] and afTect human pulmonary 
health [Pope and Dockery, 1996]. Because of its ability to 
absorb light [Horvath and Trier, 1993], soot is also recog- 
nized as an important player in the aerosol radiative forcing 
of climate at global, regional, and local scales [Menon et al., 
2002; Chung and Seinfeld, 2005; Roeckner et al., 2006]. The 
source of soot particles is the incomplete combustion of car- 
bon containing material, which means that except for natu- 
ral biomass burning all sources of soot are of anthropogenic 
origin [Penner, 1995]. The dominant removal process is wet 
deposition [Ducret and Cachier, 1992]. Soot particles can 
be transported over long distances reaching remote regions 
such as the Arctic [Clarke and Noone, 1985; Hansen and 
Nazarenko, 2004]. 

The initial composition of soot particles consists of black 
carbon and organic carbon. The precise mixture depends 
heavily on the source [Medalia and Rivin, 1982; Andreae 
and Gelencser, 2006]. While freshly emitted soot particles 
are rather hydrophobic and present in an external mixture, 
their hygroscopic qualities can change due to coagulation 
with soluble aerosols, condensation of secondary organic and 
inorganic species, and photochemical processes [Weingart- 
ner et al., 1997]. These processes arc usually referred to 
as "aging," and they determine the particle growth in re- 
sponse to ambient relative humidity and the ability to be 
activated as cloud condensation nuclei. The aging processes 
also have a profound effect on the aerosol optical properties. 
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For example, internally mixed soot shows greater absorptiv- 
ity compared to externally mixed soot. This effect on radia- 
tive properties has been studied by a number of investiga- 
tors, e.g. Chylek et al. [1995]; Jacobson [2001]; Rierner et al. 
[2003]; Schnaiter et al. [2005]; Bond et al. [2006]. Measure- 
ments show that atmospheric soot particles are internally 
mixed with other aerosol species in varying proportions, and 
that the hydrophobic portion of the aerosols decreases signif- 
icantly as the distance from the sources increases [Andreae 
et al., 1986; Levin et al., 1996; Okada and Hitzenberger, 
2001; Johnson et al., 2005; Cubison et al, 2008]. 

Since it is well recognized that soot particles contribute 
to both the direct and indirect/semi-dircct climate effect 
[Lesms et al, 2002; Jacobson, 2000, 2002b; Nenes et al, 
2002], an adequate representation of soot and its mixing 
state is sought for use in both global and regional mod- 
els, and the parameterization of soot aging is key to de- 
termining its atmospheric abundance. Many global models 
have simulated both (fresh) hydrophobic soot and (aged) 
hydrophilic soot, which can be considered as a minimal rep- 
resentation of the soot mixing state. Several of the models 
have assumed that the conversion from hydrophobic to hy- 
drophilic soot can be treated as an exponential decay pro- 
cess, with a half-life of approximately 24 h [Cooke et al, 
1999; Lohmann et al, 1999; Koch, 2001; Chung and Sein- 
feld, 2002]. This approach is a substantial simplification 
since the conversion rate depends on many dilferent envi- 
ronmental conditions. This has led to more mechanistic 
approaches, where processes such as condensation of sul- 
fate on soot particles, chemical oxidation and/or coagula- 
tion between different particle classes are explicitly modeled 
to some extent [Wilson et al, 2001; Stier et al, 2005; Tsi- 
garidis and Kanakidou, 2003]. Koch [2001] and Croft et al. 
[2005] compared different aging parameterizations in global 
models and concluded that the model results critically de- 
pend on the respective formulation. 

To better understand the soot aging process it is desir- 
able to have models that are capable of representing the 
aerosol mixing state. From a computational standpoint, if 
the aerosol mixing state can be defined in terms of A classes 
of chemical components (e.g., A = 8 with sulfate, nitrate. 
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ammonium, sea salt, hydrophobic organics, soluble organ- 
ics, black carbon, and mineral dust classes), then the mix- 
ing state is an ^-dimensional space and the size-resolved 

particle composition distribution is a multivariate function. 

Most existing aerosol models, however, represent the par- 
ticle population only as a bulk, or as a univariate function 
of a single variable, typically total mass, diameter, or simi- 
lar. To do this it is generally assumed that the population 
is either fully externally mixed, or is internally mixed with 
all particles in the same size bin or mode having the same 
mixing state. Within this framework the standard meth- 
ods are sectional, modal, and moment models. Sectional 
models (e.g. Wexler et al. [1994]; Jacobson [1997]; Adams 
et al. [1999]; Zaveri et al. [2008]) place a grid on the de- 
pendent variable space and store the number distribution 
or mass distribution (or both) on each grid cell. Modal 
models (e.g. Whitby et al. [1991]; Whitby and McMurray 
[1997]; Wilson et al. [2001]; Stier et al. [2005]; Bmkowski and 
Shankar [1995]) represent the particle distribution as a sum 
of modes, each having a log-normal (or similar) size distri- 
bution described by a small number of parameters (typically 
number, mass, and width). Moment models (e.g. McGraw 
[1997]) do not explicitly resolve the distribution, but rather 
track a few low-order moments of it. 

It is possible to extend the standard aerosol models to 
handle multivariate distributions, for example a two dimen- 
sional distribution that is a function of two species, or a func- 
tion of volume and area. Such extensions have been inves- 
tigated for sectional models [Fassi-Fihri et al., 1997], modal 
models [Brock et al, 1988], and moment models [Yoon and 
McGraw, 2004a, b]. All such models, however, require stor- 
age and computation that scale exponentially in the number 
of independent variables A. For the model we develop here 
with ^ = 20 species, fully-resolved multivariate sectional, 
modal, or moment models are infeasibly expensive. For ex- 
ample, a sectional model normally uses on the order of 8-20 
size bins to adequately resolve a univariate aerosol distri- 
bution, and even then will suffer from numerical diffusion 
[Dhaniyala and Wexler, 1996; Wu and Biswas, 1998]. An 
^-dimensional distribution would thus require 8^-20^ bins, 
which is infeasible unless A is much smaller than our 20 
species. In contrast, the particle-resolved methods devel- 
oped in this paper scale with the rmmber of particles, not 
the dimension of the space they are in. 

While traditional univariate aerosol models are too ex- 
pensive if extended to resolve multivariate aerosol mixing 
states with more than a few dimensions, there have been a 
number of extensions proposed to resolve the mixing state 
to some extent. One example of methods that somewhat re- 
solve the mixing state are the so-called source oriented mod- 
els developed by Eldering and Cass [1996] Kleeman et al. 
[1997], and Kleeman and Cass [1998] for regional scale mod- 
eling. In these models, the particles of different sources 
remain externally mixed and a number of individual size 
distributions (usually about ten) are tracked, while their 
mixing states change due to condensation of secondary sub- 
stances. However, because the main focus of their stud- 
ies was the prediction of particle mass size distributions, 
the changes in number concentrations and particle mixing 
states due to self- and hetero-coagulation of particles from 
different sources was ignored. Coagulation between aerosol 
particles is important if one is interested in predicting the 
number size distribution, especially under polluted condi- 
tions or if long residence times are considered [Zhang and 
Wexler, 2002]. Nevertheless, the source-oriented approach 
allows the attribution of pollutants to specific sources and 
is useful for designing emission control strategies [Kleeman 
and Cass, 1998; Kleeman et al., 1999]. It was used in the 
framework of a Lagrangian trajectory model, compared to 
measurements by Bhave et al. [2002] , and has been extended 
to a 3D Eulerian model [Kleeman et al, 2001; Ying et al.. 



2004; Ying and Kleeman, 2006]. The Lagrangian version de- 
scribed in Kleeman and Cass [1998] treated the mixing state 
to some extent, with fresh emissions introduced as new size 
distributions at every hour along the trajectory. 

Another variant of mixing state modeling was presented 
by Jacobson [2002a], where a total of 18 interacting aerosol 
size distributions were considered. His approach was not 
source oriented to the degree of Kleeman and coworkers, 
since anthropogenic emissions from specific sectors were not 
resolved, but primary mineral dust, searsalt, organic matter, 
and black carbon were treated. Three distributions repre- 
sented black carbon at different degrees of internal mixing. 
Coagulation between different particle classes was included, 
and 11 of the 18 particle classes were used to represent the 
mixed particles that arise due to coagulation interaction 
of two primary species. Interactions that would result in 
the formation of a particle with three diff'erent constituents 
resulted in a "mixed" particle and were not tracked fur- 
ther. Despite this considerable complexity the limitation 
remained that particles for a certain particle class and size 
were considered to be internally mixed, and the emissions 
into the primary particle categories were instantly aged. 

Riemer et al. [2003] presented an approach for mixing 
state modeling of soot using a mesoscale modal modeling 
framework. Five modes described the composition and size 
distribution of sub-micron particles, consisting of one exter- 
nally mixed soot mode, two internally mixed soot-free modes 
(containing inorganic and organic species), and two inter- 
nally mixed soot-containing modes. The last two modes thus 
represented aged soot particles, and aging occurred either by 
coagulation between modes or by condensation of secondary 
substances. While this treatment allowed the distinction be- 
tween fresh and "aged" soot, the simplifying assumption was 
made that each mode itself is internally mixed. 

Here we present a particle-resolved model, PartMC, that 
explicitly stores the composition of many individual aerosol 
particles (about 10 ') within a well-mixed computational vol- 
ume. Relative particle positions within this computational 
volume are not tracked, but rather the coagulation pro- 
cess is simulated stochastically by assuming that coagulation 
events are Poisson distributed with a Brownian kernel. 

Applying such a Monte Carlo approach for simulating the 
evolution of particle distributions dates back to Gillespie 
[1975], who developed the exact Stochastic Simulation Al- 
gorithm (see also Gillespie [1976], Gillespie [1977| and Gille- 
spie [1992]) to treat the stochastic collision-coalescence pro- 
cess in clouds. Variants of Gillespie's algorithm are widely 
used in different fields, including simulations of gene regu- 
latory networks [El Samad et al., 2005], chemical kinetics 
[Gillespie, 2007], and sintering in flames [M^e/Zs et al., 2006]. 

Since Gillespie [1975] particle-resolved methods have 
been used to study aerosols by many authors. We do 
not attempt to give a comprehensive literature survey 
here. Babovsky [1999] and Eibeck and Wagner [2001] 
developed the Mass Flow Algorithm with variable com- 
putational/physical particle ratios, Kolodko and Sabelfeld 
[2003] gave relevant error estimates, and Debry et al. [2003] 
coupled it to evaporation and condensation. Somewhat 
similarly, Laurenzi et al. [2002] and Alfonso et al. [2008] 
(based on ideas from Spouge [1985]) stored the number 
of particles with identical composition to reduce memory 
usage and computational expense while using Gillespie's 
method. Guias [1997] studied convergence of stochastic 
coagulation to the Smolukowski equation. Efendiev and 
Zachariah [2002] investigated enclosures within aerosols us- 
ing a particle-based method, while Maisels et al. [2004] used 
particle methods with simultaneous nucleation, coagulation, 
and surface growth. 

While not focused on aerosol simulations, much re- 
cent work has investigated efficient simulation methods for 



Preprint generated on Thursday 4* September, 2008 — Submitted to Journal of Geophysical Research 



RIEMER ET AL.: PARTICLE-RESOLVED MIXING STATE MODELING 



3 



reaction-type Markov processes. Gillespie [2001] developed 
the tau-leaping method for efficient generation of many 
events with near-constant rates, with extensions by Gillespie 
and Petzold [2003] ; Rathinam et al. [2003]; Cao et al. [2006] 
and others, including for multiscale systems with scale sepa- 
ration in Cao et al. [2005]. Multiscale variants of Gillespie's 
Stochastic Simulation Algorithm have also been developed 
by E et al. [2007]. Gibson and Brack [2000] developed the 
Next Reaction Method for efficient exact sampling, which 
stores and reuses event calculations for efficiency. Ander- 
son [2007] and Anderson [2008] developed efficient simula- 
tion algorithms based on the Next Reaction Method and the 
tau-leaping method. 

For the large number of particles in the simulations in this 
paper, we used an efficient approximate coagulation method, 
as described in Section 4 below. This used a binned sam- 
pling method to efficiently sample from the highly multiscale 
coagulation kernel in the presence of a very non-uniform 
particle size distribution, implemented with a multi-event- 
per-timestep sampling of the coagulation events. Multi-rate 
versions of Gillespie's method have been developed previ- 
ously by Cao et al. [2005] and E et al. [2007], but relied 
on scale separation to average slow event rates over fast 
timcscalcs. The method used here does not accelerate rare 
events but it does accelerate the generation of events without 
scale separation, as needed for the smoothly varying coagu- 
lation kernels and particle size distributions. The PartMC 
coagulation method has storage cost proportional to the 
number of physical particles, computational cost for evapo- 
ration/condensation proportional to the number of particles, 
and computational cost for coagulation proportional to the 
number of coagulation events. 

PartMC was coupled with the new state-of-the-art aerosol 
chemistry model MOSAIC [Zaveri et al, 2008], which simu- 
lates the gas- and particle-phase chemistries, particle-phase 
thermodynamics, and dynamic gas-particle mass transfer in 
a deterministic maimer. MOSAIC treats all the important 
aerosol species, including sulfate, nitrate, chloride, carbon- 
ate, ammonium, sodium, calcium, primary organic mass 
(POM), secondary organic mass, black carbon (BC), and 
inert organic mass. The coupled model system, PartMC- 
MOSAIC, accurately predicts number, mass, and composi- 
tion size distributions, and is therefore suited for applica- 
tions where any or all of these quantities are required. 

Simulating all particles explicitly in a population of 
aerosol completely eliminates any errors associated with nu- 
merical diffusion. As a result, the highly accurate treat- 
ment of aerosol dynamics and chemistry makes PartMC- 
MOSAIC suitable for use as a numerical benchmark of 
mixing state for more approximate models. It can also 
be applied to different environments going beyond the ex- 
ample of clear-sky photochemistry shown in this paper, 
including the in-cloud processing of aerosol, and it can 
be used to accurately estimate quantities that depend on 
the mixing state, such as cloud condensation rmclci spec- 
tra and optical properties, which we will address in a 
forthcoming paper. The current version of PartMC is 
available under the GNU General Public License (GPL) 
at http:/ /www. mechse.uiuc.edu/research/mwest/partmc/, 
while the MOSAIC code is available upon request from 
R. A. Zaveri. 

This manuscript is organized as follows. In Section 2 we 
write the governing equations for the coupled gas-aerosol 
box model and discuss the approximations needed by this 
model of the physical system. The numerical approximation 

to the governing equations is given in Section 3, where we 
introduce the particle-resolved aerosol model PartMC and 
describe how it is coupled to the gas- and aerosol-chemistry 
code MOSAIC. In Section 4 we give the efficient coagulation 
algorithm used by PartMC and verify its performance nu- 
merically. Finally, Section 5 demonstrates the capabilities 
of this new model approach by focusing on the evolution 



of the mixing state of soot particles in an idealized urban 
plume scenario. The main contributions of this paper are: 1) 
an accelerated stochastic coagulation method for multiscale 

kernels, 2) the coupling of a particle resolved model with a 
gas- and aerosol-chemistry code, and 3) an initial study of 
the soot mixing states present in a typical polluted urban 
environment. 

2. Coupled aerosol-gas governing equations 

We consider a Lagrangian parcel framework where we 
simulate the evolution of aerosol particles and trace gases 
in single parcel (or volume) of air moving along a specified 
trajectory. In addition to coagulation and aerosol and gas 
chemistry, the model treats prescribed emissions of aerosols 
and gases, and mixing of the parcel with "background air". 

The evolution of atmospheric aerosols is extremely com- 
plex, involving interactions between fluid-transport and 
micro-scale properties of the aerosol particles during coag- 
ulation and condensation. A fully resolved simulation of 
fluid-aerosol interaction cannot capture a large enough sys- 
tem to determine macro-scale aerosol distribution proper- 
ties. Instead of fully resolved models, it is usual to use a 
box model locally and store the size distribution of aerosol 
in a certain physical volume without storing the positions of 
the particles in throe dimensions. Particle interactions such 
as coagulation arc then represented stochastically by a ker- 
nel K that defines coagulation probability rates for particles 
depending on their sizes and compositions. This representa- 
tion assumes that the mixing timescale is significantly faster 
than the coagulation and condensation timescales. Equiv- 
alently, we are assuming that the aerosol processes are al- 
most Markovian and so a memoryless model is a good ap- 
proximation to the physics. This assumption underlies all 
sectional and modal models in use today, as well as our 
particle-resolved model PartMC, and means that we need 
only consider the aerosol- and gas-species densities. 

An aerosol particle contains mass /Xa > of species a, for 
a = 1, . . . ,A, so that the particle composition is described 
by the j4-dimensional vector p, € R'*. jUaii (m^) is the total 
wet mass of the particle, and fidiy = Mali — A'H20 (m^) is 
the total dry mass. The cumulative aerosol number distri- 
bution at time t and constituent masses fl G is N{il, t) 
(m~^), which is defined to be the number density of aerosol 
particles that contain less than /ia mass of species a, for all 
a = 1, . . . ,A. The aerosol number distribution at time t and 
constituent masses ft € R'* is n{jl,t) (m~^kg~'^), which is 
defined by 

nip,t)= ' . 1 

dij,idij,2 . . . OUA 

The concentration of trace gas phase species i at time t 



is given by gi{t) (molm ), for 



. , G, so the trace 



gas phase species concentrations are the G-dimensional vec- 
tor g{t) G R'-'. We assume that the aerosol and gas species 
are numbered so that the first C species of each undergo 
gas to particle conversion, and that they are in the same 
order so that gas species i converts to aerosol species i, for 
i = l,...,C. 

The environment is described by temperature T{t) (K), 
pressure p(t) (Pa), relative humidity RH(t) (dimensionless) , 
and dry density Pdry(<) (kgm^^). For the simulation in Sec- 
tion 5 the air temperature is prescribed as a function of time, 
while the air pressure and water mixing ratio are kept con- 
stant and the relative humidity and dry density are updated 
accordingly. 

We assume that we are modeling a vertical slice of a well- 
mixed boundary layer during the day and a slice of the resid- 
ual layer during the night, always surrounded to the sides 
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and above by background air. The height of the boundary 
layer is given by H{t) (m). We denote by Adii,horiz(i) (s~^) 
the horizontal dilution rate with the prescribed background 
gas and aerosol, and by Adii,vert(^) (s~^) the vertical dilu- 
tion rate that represents entrainment of a growing boundary 
layer. The total dilution rate A|^(^) (s~^) is then given by 



Adil(t) = Adil,horiz(t) + Adil,vert(^) (2) 

\ ' H{t) dt 



Adil,vert = /entrain (t) max( 0, '^^}f^ \ : (3) 



where vertical entrainment only occurs for increasing bound- 
ary layer heights. The indicator /entrain (i) is 1 when the 
modeled air parcel is within the boundary layer and so en- 
trainment is possible, and is when the air parcel is in the 
residual layer. 

To obtain the mean evolution in the large-number limit 
we neglect correlations between the number of particles of 
different sizes [Gillespie, 1972] and thus obtain the classi- 
cal Smoluchowski coagulation equation [von Smoluchowski, 
1916a, b], which for a multidimensional aerosol distribution 
with gas coupling is 

dn{ii,t) 1 ^ ^. 

X n(/Lt', t) n{fl — ,t) d/i'i d/Lt2 ■ ■ • ^Ma 

(4a) 

roc roc roo 

- / •■■/ K{il,fl') 
JO Jo Jo 

X n{iJ,t)n{ff ,t) dfjfidf/2 ■ ■ ■ (4b) 

+ ncmit(P,t) (4c) 

+ Adii {t) (jihi^ck ifi,t) - n{jl, t)J (4d) 



(4e) 
(4f) 
(4g) 



~ ^ ( Ci/i(M, 9, t) n{fl, t) 



d 



dfic+i 
1 



Pdry(t) dt 



Cw/w(/i, g, f)n(/I, t) 
n{U, t) 



iry(i)^ 



dgijt) 
dt 



■■ 5emit,i it) + Adil {t) {gha.ck,i (*) - 9i (*)) (4h) 

+ Ri{a) (4i) 

/•CO noo roo 

- / ••• / ii{fi,9,t) 

Jo Jo Jo 

X n{fl, t) d/ui d/i2 ■ ■ ■ dfj,A (4j) 

1 dpdry (i) 



Pdry(i) at 



(4k) 



The integro-differential equation (4) must be augmented 
with appropriate boundary conditions. This is done on 
physical grounds to ensure that the constituent masses 
of particles cannot become negative. In equation (4), 
K{jli,p,2) (m^s~^) is the coagulation rate between par- 
ticles with constituent masses jl\ and /i2, riemit (/*,*) 



(m-^kg-^s- 



is the number distribution rate of aerosol 



emissions, nbacii(A', ^) (ni~ kg~ ) is the background num- 
ber distribution, Ii{p,g,t) (mols^^) is the condensation or 
evaporation flux of gas species i (with g, t) the flux for 

water), Ci {m^ moP ) is the conversion factor from moles of 
gas species i to mass of aerosol species i (with Cw the fac- 
tor for water), 3emit,i(i) (molm^ s^^) is the emission rate 
of gas species i, ghs.c'k.iit) (molm^"') is the background con- 
centration of gas species i, and Ri{g) (molm^ s^^) is the 
concentration growth rate of gas species i due to gas chem- 
ical reactions. Many of the rates, coefficients and functions 



also depend on the environmental conditions, but we have 
not written this dependence explicitly. 

3. Particle-resolved aerosol models 

3.1. PartMC aerosol state representation 

We consider a Lagrangian parcel with volume V (m*), 
also called the computational volume. We represent the 
aerosol state by storing Nmg particles in this volume, writ- 
ten n = {p} ,fp , . . . ,fl'^^'^), where the particle order is 
not significant. Each particle is an ^-dimensional vec- 
tor fT G R"* with components (/il, //2, • • • , Ma); so /i^ is 
the mass of species a in particle i, for a — 1,. . . ,A and 
i — 1, . . . , A^MC- In the notation of Debry et al. [2003] for 
the Mass Flow Algorithm, we are taking {iUi/yi){t) = 1, 
which means one computational particle per physical parti- 
cle. While we track every particle within the computational 
volume V, we regard this volume as being representative of 
a much larger air parcel. For example, in Section 5 we use 
a computational volume on the order of a few cubic cen- 
timeters but take this to be approximating the state of the 
well-mixed boundary layer during the day and the residual 
layer during the night. 

The simulation of the aerosol state proceeds by two mech- 
anisms. First, the composition of each particle can change, 
changing the components of the vector fT as species con- 
dense from the gas phase and evaporate to it, for example. 
Second, the population II can have particles added and re- 
moved, either by emissions, dilution or coagulation events 
between particles. 

The representation of the aerosol as a finite collection of 
particles 11 in a volume V is very flexible, as other properties 
can easily be stored for each particle, such as fractal dimen- 
sion, electric charge, age since emission, etc. In the present 
paper wo store the number of coagulation events undergone 
by each particle to produce Figure 15. 

3.2. PartMC emissions 

Because we are using a finite number of particles to ap- 
proximate the current aerosol population, we need to add 
a finite number of emitted particles to the volume at each 
timestep. Over time these finite particle samplings should 
approximate the continuum emission distribution, so the 
samplings at each timestep must be different. We assume 
that emissions are memoryless, so that emission of each par- 
ticle is uncorrelated with emission of any other particle. Un- 
der this assumption the appropriate statistics are Poisson 
distributed, whereby the distribution of finite particles is 
parametrized by the mean emission rate and distribution. 

Consider a number distribution production rate riemit (m, t) 
(m~^ kg~^ s~^), a volume V (m^), and a timestep At (s). 
The omissions over the timestep from time to to ti = to + At 
are given by 



«emit iP;to,tl) = / Uemit (M, t) dt 
J to 



{tl — to)«'emit(At, to) 



(5) 
(6) 



for which we use the first-order approximation above. To 
obtain a finite Poisson sample of the distribution n(/2) — 
ncinit{fl',to,ti) (m^'^kg^'^) in the computational volume V 
we first see that the mean number N{n, V) of sampled par- 
ticles will be 



iVmean(n,V^) : 



oc poo 



n{fT)V djjLi . . . df^A- 



(7) 

The actual number S of emitted particles added in a 
timestep will be Poisson distributed, written S ~ Pois(A), 
for mean A = A^mean(n, V), so that 



Prob(S' = k) 



k\ 



for fc e ! 



(8) 
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A Poisson sampling Ilsamp of the number distribution n{jT) 
in volume V, written Ilsamp ~ Poisdist(w, 1^), is a finite se- 
quence of particles given by 



n, 



samp - {jl ,fl ) (9a) 

I, V)) (9b 

fors = l,...,S', (9c) 



S'~P0is(iVmean(n,l/)) (9b) 

n{il)V 



Armean(n, V) 



where (9c) means that each particle has a composition drawn 
from the distribution specified by n{jT). 

3.3. PartMC dilution 

As with emissions, we must also obtain a finite sampling 
of background particles that have diluted into our computa- 
tional volume during each timcstcp. In addition, some of the 
particles in our current sample will dilute out of our volume 
and will be lost, so this rrmst be sampled as well. We assume 
that dilution is memorylcss, so that dilution of each parti- 
cle is uncorrelated with the dilution of any other particle 
or itself at other times, and that once a particle dilutes out 
it is no more likely to re-enter than any other background 
particle. 

Let the background particle distribution be nback(/i, t) 
(m~^kg~^), the computational volume be V (m^), and the 
timestep be At (s). The distribution of particles that dilute 
from the background into the volume V between times to 
and ti =to + At is ndii(At; to,ti), where ndii(At; to, t) satisfies 

dndn{P^to,t) ^ ^nback(/i, t) - ndii(M; to, t)^ (10a) 

ndii{i2;to,to) = (10b) 



n is the sequence of particle compositions 
V is the computational volume 
g is the gas concentrations 
t = 

while t < tfinai do: 
t = t + At 

update temperature T{t), pressure p{t), relative 
humidity RH{t), dry density pdry(t), and 
mixing height H{t) 

V{t)=V{t-AtY-^^^^ 

add Atpemit(i) + At Adii(t)^ffback(^) - g{t)J to g 

randomly choose iVioss ~ Pois(At AdiiATMc) and 
remove Mogg randomly chosen particles from 

n 



add a sample of Poisdist i^^du At nback(-, t), Vj to n 
add a sample of PoiSdist ( At Tiomit {■,t),V] to 11 



perform one At-timestep of coagulation for n with 
the PartMC algorithm in Figure 2 

integrate the system of coupled ODEs (4e) 
and (4i)-(4j) with MOSAIC for time At 

end while 



We use the first-order approximation given by 

ridiiiP; to, ti) ~ (ti - to) Adii(to) nback(/t, to) (11) 

A discrete sampling of ndii(/i; to, ti) is then given by Ildii ~ 
Poisdist(ndii(/t), V), as in (9). 

If wc start the timestep at time to with the particle pop- 
ulation n, then each particle in 11 has probability p(to,ti) 
to be lost by dilution during the timestep, where p(to,t) 
satisfies 



dp{to,t) 

dt 

pita, to) = 



Adii(t)(l-p(to,t)) (12a) 

(12b) 



We use the first-order approximation given by 
p{to,ti) ^ (ti -to)Adii(to) 



(13) 



We denote the binomial distribution for number n and prob- 
ability p by B{n,p). The number of particles lost from II 
between times to and ti = to -|- At is then given by Moss 
which is distributed as 



Moss - B(^iVMC,p(io,tl)) 

w Pois(MMcp(to,ti)j 



(14) 
(15) 



We approximate the binomial distribution with a Poisson 
distribution as above, which converges as At — > for fixed 

Nmc- As each particle has equal probability to be lost due 
to dilution, we can sample Mioss and then choose Moss par- 
ticles uniformly from 11 to be removed. 

3.4. Coupled PartMC-MOSAIC method 

We couple the stochastic PartMC particle-resolved 
aerosol model to the deterministic MOSAIC gas- and 
aerosol-chemistry code in a time-splitting fashion to obtain a 
complete discretization of the governing equations (4). The 
aerosol distribution n(/i, t) is represented by Nmc particles 
in a computational volume V, as described above, while 
the gas vector g{t) stores the gas concentrations. Equa- 
tions (4a)-(4d) are solved stochastically by the PartMC 
code, and the gas equations (4h)-(4j) together with the cou- 
pling term (4e) are integrated deterministically by the MO- 
SAIC code. The terms (4g) and (4k) are implemented deter- 
ministically by updating V and g by the density change or its 
inverse, as appropriate. The full coupled PartMC-MOSAIC 
algorithm is given in Figure 1. 

The current version of MOSAIC treats all the locally 
and globally important aerosol species including SO4, NO3, 
CI, CO3, MSA (methanesulfonic acid), NH4, Na, Ca, other 
inorganic mass, BC, and POM, and secondary organic 
mass. It consists of four computationally efficient mod- 
ules: 1) the gas-phase photochemical mechanism CBM-Z 
[Zaveri and Peters, 1999]; 2) the Multicomponcnt Taylor 
Expansion Method (MTEM) for estimating activity coeffi- 
cients of electrolytes and ions in aqueous solutions [Zaveri 
et al, 2005b]; 3) the Multicomponcnt Equilibrium Solver for 
Aerosols (MESA) for intra-particle solid-liquid partitioning 
[Zaveri et al., 2005a]; and 4) the Adaptive Step Time-split 
Euler Method (ASTEM) for dynamic gas-particle partition- 
ing over size- and composition-resolved aerosol [Zaveri et al., 
2008]. The version of MOSAIC box-model implemented 
here also includes a treatment for secondary organic aerosol 
(SOA) based on the SORGAM scheme [Schell et al, 2001]. 



Figure 1. Coupled PartMC-MOSAIC algorithm. 
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4. PartMC coagulation algorithm 
4.1. Stochastic coagulation simulation 

If wc have Nmc particles then there are iVMc(AfMC — l)/2 
possible coagulation events, with the probability rate of a co- 
agulation between particles i and j in a volume V given by 
K{fr , fl^) /V for the coagulation kernel K(jl\p-') (rn'^s^^). 
The only difficulty with stochastic coagulation is generating 
a sequence of coagulation events, each consisting of a pair of 
particles (i, j) that coagulate and a time At until the coagu- 
lation occurs. Once a coagulation event is determined then 
we simply remove the particles i and j from the population 
n, add a new particle with composition = fT -\-pP , and 
advance the time by At. 

The standard stochastic simulation algorithm for this sys- 
tem is due to Gillespie [1975] and is based on the observa- 
tion that the probability density for the time until the next 
coagulation event is 



P(At) 



_tot -iftotAt/V 
V 



(16) 



where JsTtot = l^i<,' ^ifi^^i^^) is the total rate. We can thus 
generate an elapsed time by sampling the probability density 
function (16). The conditional probability that the coagu- 
lation event that occurred was between particles i and j is 
then 

P(..|At)^^, (17) 

and this can be sampled to determine which particles coag- 
ulated, and then the coagulation event can be performed. 

Gillespie's method has the advantage that it generates 
exact realizations of the stochastic coagulation process. It 
faces two main difficulties in practice, however. First, the 
total rate Ktot continually changes as coagulation events oc- 
cur and particle compositions change due to condensation. 
Computing a reasonable estimate of this parameter quickly 
becomes exceedingly expensive, and approximations made 
to speed up this estimate introduce errors that arc difficult 
to estimate and control. Second, while sampling (16) is very 
cheap, sampling (17) can be expensive for complex kernels. 
The two main methods arc use of the cumulative distribu- 
tion function, which scales badly in the number of particles 
and is thus too expensive for large particle rmmbers, and 
use of acccpt-rcjcct. While accept-reject scales well as the 
number of particles grows, it is very inefficient if the ker- 
nel Kij is highly non-uniform, as is unfortunately the case 
for many physically relevant aerosol distributions. Despite 
these difficulties, Gillespie's method is by far the most com- 
monly used method in practice, with many slight variants 
appearing in the literature (for example, see Efendiev and 
Zachanah [2002]; Krms et al. [2000]; Garcia et al. [1987]; 
Fichthorn and Weinberg [1991]). 

To avoid these two difficulties wc fornmlatc an improved 
method. Wc use a fixed timestep method and wc develop a 
binned acceptance procedure. The use of a fixed timestep 
removes the need to know iftot, albeit with the introduc- 
tion of some error. This fixed timestep also makes it easy 
to integrate the coagulation with other physics and chem- 
istry using a time-splitting scheme. The binned sampling 
method means that we are not subject to slow-downs from 
non-uniform kernels. 

4.2. Fixed-timestep stochastic coagulation 

Wc choose a fixed timestep At and in each timestep 
choose Attest particle pairs to test. We then generate Attest 
random particle pairs uniformly and for each pair (i, j) we 
accept a coagulation event with probability 



-P(«>i) = 1 - exp 



-A'(/r,/2^)At7VMc(iVMC-l) 



V iVtest 
ii:(/I\/I^)AtiVMc(iVMC-l) 



V 



(18) 



In the limit At — » this generates an exact realization of 
the stochastic coagulation process, and for finite At intro- 
duces a discretization error. The number Attest should be 

chosen large enough that P{i,j) < 1 for all pairs i,j and 
for convergence it must remain bounded away from zero as 
At 0. This is similar to the sampling technique used in 
Debry et al. [2003]. 

The efficiency of the method, as with any procedure of 
accept-reject type, is greatest when the maximum value of 
P{i,j) is as close as possible to 1. To ensure this we choose 



A/'test = ["A'maxAtA/-Mc(AfMC - l)/vj , 



(19) 



where JCmax = maxj,j K^fT ,fp) is the maximum kernel value 
and \x\ is the least integer greater than x. In practice 
we take /Cmax to be a cheaply computable upper-bound 
for K{fr,pP), which slightly increases the accuracy of the 
method and is much cheaper. 

The fixed timestep method thus cleanly resolves the dif- 
ficulties with Gillespie's method to do with the need to de- 
termine /Ctot. It still has the problem, however, that if the 
kernel K{ir,pP) is very non-uniform then the acceptance 
procedure will be very inefficient. To fix this, we adopt a 
binned approach. 

4.3. Binned stocheistic coagulation 

For coagulation kernels of physical interest, such as those 
arising from Brownian motion or gravitational settling, the 
kernel K{p,p') is highly multiscale, with many orders of 
magnitude difference between the highest and lowest rates. 
This is a problem for the sampling procedure outlined in the 
previous section, because A/test will be very large and so we 
will have to reject many events using (18) for each accepted 
event. 

To accelerate this procedure we take advantage of the 
fact that the kernel K{p',il^) is not random in its non- 
uniformity, but rather depends primarily on the diameter 
of the particles. This means that if pairs i,j and k,£ are 



divide diameter axis into bins as for a sectional model 
A'Mc(b) is the number of particles in bin b 
fl{b, i) is the mass vector of the i-th particle in bin h 
^max (61,62) is a precomputed upper bound on the kernel 

for any particles from bins 61 and 62 
At is the timestep 
for all bin pairs (61,62) do: 

ATevent = ATmC (61 )ArMC (62)72 
Attest = |"/fmax(6l,62) AtA/"event/V^] 

for A/test repetitions do: 

randomly choose particles ii and 12 uni- 
formly in bins 61 and 62 

if 12 = K[fl{b\,i\),jl{h2,i2)) 
randomly choose r uniformly in [0, 1] 
if r < Ki2 AtNevem/iNtestV) then: 

coagulate the two particles, updating 
the arrays N{b) and jl{b, i) 

end if 



end for 



end for 



Figure 2. PartMC coagulation algorithm. 
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similar, so that the diameters of particles i and k are close, 
as are the diameters of particles j and £, then K(p^,fp) w 
K(fl'',fl'). We thus group particles into bins sorted by di- 
ameter and wc use the acceptance procedure (18) for each 
pair of bins separately. This binned approach ensures that 
all particle pairs under consideration in a particular itera- 
tion have similar coagulation rates, and hence the procedure 
will have a high proportion of acceptances. Use of a binned 
version of the fixed timcstcp algorithm means that the num- 
ber of samples (19) done per pair of bins is automatically 
adapted to the number of particles in those bins. It also 
allows us to pre-compute the ifmax values for each bin pair. 
The resulting algorithm is shown in Figure 2. 

The primary disadvantage of using a binned sampling 
procedure is in code complexity, as the bin structures of 
particles with similar sizes need to be constructed and main- 
tained. This also adds a small amount of computational 
overhead to the coagulation routine, which is far outweighed 
by the enormous efficiency gains. We should note that the 
binned sampling procedure introduces no error in the sim- 
ulation and is a pure efficiency gain. For typical aerosol 
profiles the binned procedure gives about two to four orders- 
of-magnitude speedup in computational time, as quantified 
in Section 4.4. 

As a result of continuous coagulation the total number of 
simulated particles decreases with time. To maintain sufli- 
cient resolution and adequate statistics we double the num- 
ber of particles whenever the number of particles becomes 
less than half of the original particle number. This corre- 
sponds to a doubling of the computational volume. When 
we include emissions of particles, the particle number even- 
tually becomes too large and we run into computational lim- 
its. To avoid this we halve the number of particles when we 
have more than twice the original number of particles, which 
corresponds to a halving of the computational volume. 

For some kernels, such as the Brownian kernel used in 
Section 5, the kernel is primarily dependent on the particle 



2- 10" 




dry diameter (/jm) 

Figure 3. Comparison of the stochastic particle-resolved 
method using 10® particles (circles) against a sectional so- 
lution (lines) to the Smoluchowski equation for the Brow- 
nian kernel according to Jacobson [1999] . 



1.5-10' 




elapsed time (hours) 

Figure 4. Comparison of the stochastic particle-resolved 
method using 10"" particles (circles) against the analyt- 
ical solution (linos) for a simulation with only constant 
mean rate emissions and dilution. 

diameters but also depends on particle density. We could 
store the particles sorted into a 2D array per-diameter-bin 
and per-density-bin, but the density variation is bounded 
and small enough that it is still reasonably efficient to store 
them only per-diameter-bin and to compute -fCmax to take 
the maximum density variation into account. 

To enable efficient coagulation, the particle array If is 
stored as an array of pointers to partially-filled particle ar- 
rays, one per diameter-bin. Insertions into bin arrays are 
performed at the end of the currently filled area and dele- 
tions from the middle are followed by a shift of the last 
element into the gap, ensuring full packing of each bin ar- 
ray at all times. Each diameter-bin array is reallocated to 
twice its existing size when necessary or half its existing size 
when possible. This gives constant-time random access at 
the cost of 0(log AA^Mc) reallocations and at most twice the 
minimal memory usage. 

4.4. Verification of the PartMC coagulation algorithm 

For verification of the PartMC stochastic coagulation 
method we compared PartMC using 10'' particles against a 
sectional solution to the Smoluchowski equation [Bott, 1998] 
for a Browruan kernel [Jacobson, 1999]. For Figure 3 wc used 
two overlapping log-normal modes as the initial condition 
and the results show that we have excellent agreement for 
the number and mass distributions for this test case, which 
is representative of the simulation in Section 5. At the very 
largest sizes there is some noise in the particle-based mass 
distribution, as each individual particle has significant mass 
at these sizes. This noise could be reduced by averaging 
several simulations in a Monte Carlo fashion, or by using 
a variable number of physical particles per computational 
particle, as in the Mass Flow Algorithm [Babovsky, 1999; 
Eibeck and Wagner, 2001]. We do not consider this noise to 
be significant enough for the study in this paper to require 
amelioration. 

For the Brownian kernel in Figure 3 the use of the binned 
stochastic coagulation algorithm of Section 4.3 improved the 
accept rate from 0.95% to 86%, requiring more than 90 times 
fewer kernel evaluations. For a more non-uniform gravita- 
tional kernel, such as found in cloud-aerosol simulations, the 
binned algorithm increased the accept rate from 0.007% to 
86%, a reduction of over 12,000 times in the number of ker- 
nel evaluations (not shown in a figure). 

Figure 4 compares the stochastic treatment of emissions 
and dilution using 10^ particles against the analytical solu- 
tion for constant mean emission and dilution rates. This test 
also shows excellent agreement. We thus see that PartMC- 
MOSAIC is performing emissions, dilution, and coagulation 
accurately, and the chemistry modeling is of similar accu- 
racy to that in current state-of-the-art sectional and modal 
aerosol models. 
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Table 1. Initial concentrations of gas phase species and average gas phase emissions. The emissions represent 
area emissions and are averaged over the 12-hour emission period. We obtain the volume emission rate Semit(*) 
in equation (4h) by dividing by the mixing height H{t). 



MOSAIC species 


Concentration in ppb 


Emissions in nmol m ^ s ^ 


NO 


0.1 


31.8 


NO2 


1.0 


1.67 


HNO3 


1.0 




O3 


50.0 




H2O2 


1.1 




CO 


21 


291.3 


SO2 


0.8 


2.51 


NH3 


0.5 


6.11 


HCl 


0.7 




CH4 


2200 


_ 


C2H6 


1.0 


_ 


HCHO 


1.2 


1.68 


CH30H 


0.12 


0.28 


CH300H 


0.5 


_ 


ALD2 


1.0 


0.68 


PAR 


2.0 


96 


AONE 


1.0 


1.23 


ETH 


0.2 


7.2 


OLET 


2.3- 10-2 


2.42 


OLEI 


3.1 • 10-4 


2.42 


TOL 


0.1 


4.04 


XYL 


0.1 


2.41 


ONIT 


0.1 




PAN 


0.8 




RCOOH 


0.2 




ROOH 


2.5 • 10^2 




ISOP 


0.5 


0.23 


ANOL 




3.45 



Table 2. Initial and emitted aerosol distribution parameters, as in equation (23). The initial aerosol distribution 
is also used as the background aerosol distribution. E is the area source strength of particle emissions. Dividing E 
by the mixing height H{t) and multiplying by a normalized composition distribution gives the number distribution 
emission rate riemit in equation (4c). 





TV in 


m 


Dgn in fim 




Composition 


Initial mode 1 


3.2 


■10« 


0.02 


1.45 


50% (NH4)2S04, 50% POM 


Initial mode 2 


2.9 


•109 


0.116 


1.65 


50% (NH4)2S04, 50% POM 




B in m ^ s 


Dgn in ixm 


<^s 


Composition 


Meat cooking 


9- 


10" 


0.086 


1.9 


POM 


Diesel vehicles 


1.6 


• 108 


0.05 


1.7 


10% POM, 90% BC 


Gasoline vehicles 


5 • 


10^ 


0.05 


1.7 


66% POM, 34% BC 



5. Application of PartMC-MOSAIC to an 
idealized urban plume scenario 

5.1. Setup of case study 

For this study we considered an idealized urban plume 
scenario. We tracked the evolution of gas phase species and 
aerosol particles in a Lagrangian air parcel that initially con- 
tains background air and is advected over and beyond a large 
urban area. The simulation started at 06:00 local standard 
time (LST), and during the advection process, primary trace 
gases and aerosol particles from different sources were emit- 
ted into the air parcel for 12 hours. After 18:00 LST, the 
emissions were switched off, and the evolution of the air 
parcel was tracked for another 12 hours. 

Initial gas-phase concentrations and emissions were 
adapted from the Southern California Air Quality Study 
(SCAQS) simulation (August 26-29, 1988 period) in Zaveri 
et al. [2008], and are listed in Table 1. Note that while gas 
phase emissions in the simulation varied with time. Table 1 
gives only the average over the emission period. The initial 
particle size distribution, which was identical to the back- 
ground aerosol distribution, was bimodal with Aitken and 
accurrmlation modes [Jaenicke, 1993]. We assumed that it 
consisted of (NH4)2S04 and primary organic mass (see Ta- 
ble 2). We considered three different types of carbonaceous 
aerosol emissions: 1) meat cooking aerosol, 2) diesel vehi- 



cle soot, and 3) gasoline vehicle soot. The parameters for 

the size distributions of these three emission categories were 
based on Eldering and Cass [1996], Kittelson et al. [2006a], 
and Kittelson et al, [2006b] , respectively. The emission rates 
and the compositions were adapted from the California Air 
Resources Board database [California Air Resources Board, 
2007]. 

For simplicity in this idealized study, the particle emis- 
sions strength and their size and composition were kept con- 
stant with time. Furthermore, we assumed that the particles 
from these sources were emitted as fully-internal mixtures of 
the species listed in Table 2, since to date the mixing state 
of particle emissions is still not well characterized. Hence, 
it is difficult to justify a more sophisticated treatment for 
the emissions. However, once this data is available it will be 
straightforward to implement the information accordingly in 
our model. Sea salt, biomass burning and mineral dust par- 
ticles as well as particles from biological sources (e.g. pollen) 
were not treated in this test case. 

Before we discuss the results on aerosol mixing state in 
detail we provide the context for the conditions in our case 
study with Figures 5 to 8. We did not attempt to sim- 
ulate a specific episode or trajectory for the Los Angeles 
basin (as was done in Kleeman et al. [1997]), but rather an 
idealized urban plume scenario, with conditions that were 
consistent with a polluted environment. The temperature, 
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Figure 5. Time series of temperature, relative humidity, and mixing height over the course of the 24-hour simulation. 



relative humidity, and mixing height along the trajectory 
were adapted from spatially-averaged values from the Los 
Angeles Air Basin simulation of Zaveri et al. [2008] and ref- 
erences therein. The temperature and mixing height were 
prescribed as functions of time, while the pressure and water 
mixing ratio were kept constant and the relative humidity 
and dry density were updated accordingly. The variation 
of these parameters is shown in Figure 5. The relative hu- 
midity started at 95%, then decreased to 53% during the 
day and increased again to 94% during the following night. 
As we show below, the diurnal cycle of the ambient condi- 
tions impacted the thermodynamic equilibria and the phase 
states of the particles. 

An increase of the mixing height during the morning 
caused dilution of the gas and aerosol concentrations within 
the air parcel and was accompanied by entrainment of back- 
ground air, as discussed in Section 3.2. We also considered 
dilution due to horizontal turbulent diffusion, using a first- 
order dilution rate of 1.5 • 10~^ms~^. 
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Figure 6. Time series of selected gas phase species. 



06:00 



We resolved the total aerosol distribution with 10^ Monte 
Carlo particles initially. The corresponding initial total 
number density was N — 6.1 • lO'' m^"' and so the com- 
putational volume was initially V = Nmc/N = 16 cm^. It 
remained between V = 8cm^ and V = 17cm^ for the du- 
ration of the run as particle number ATmc and number den- 
sity N changed due to emissions, dilution, and coagulation. 
The number of particles remained between Nmc = 199799 
and Nmc = 60655. The timestep used for this simulation 
was At = 1 minute. While better estimates of the system 
statistics could be obtained with multiple realizations, we 
found a single run to give reasonable results in this case, 
as demonstrated in Figures 3 and 4 and discussed in Sec- 
tion 4.4. Although not shown here, runs with different ran- 
dom initialization gave essentially the same results. 

To quantify the impact of coagulation we performed two 
runs, one base case including coagulation as described above. 
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Figure 7. Time series of total mass densities of selected 
aerosol species: Mnos, -^nh4, Mpom, Afso4, Mbc, and 

MsOA. 
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and one case without coagulation. Otherwise, the conditions 
for the two runs were identical. 

5.2. Gas species evolution 

Figure 6 shows the evolution of seleetcd gas phase species 
undergoing a diurnal cycle typical for a photochemistry 
episode under polluted conditions. During the daytime we 
observed a considerable production of O3, reaching a max- 
imum value of 144 ppb at 16:09 LST. The NO2 concentra- 
tion increased up to 33 ppb during the time that NOx was 
emitted, and decreased after 18:00 LST due to dilution and 
chemical reactions after the emissions had stopped. IfNOa 
reached 17 ppb and contributed to the formation of am- 
monium nitrate in the particle phase. NH3 levels reached 
6.6 ppb during the daytime and later vanished due to gas- 
to-particle conversion. HCHO was both emitted and chem- 
ically produced with a maximum value of 12 ppb at 17:59 
LST. 

5.3. Bulk aerosol evolution 

Figure 7 shows time series of the bulk aerosol concentra- 
tions. We observe a pronounced production of ammonium 
nitrate, reaching nitrate concentrations of up to 26/igm~^ 
and ammonium concentration of 10.3 /igm^'* in the late af- 
ternoon. Sulfate concentrations increased from 4.1/igm~^ 
to 6.8/Ltgm~^ due to condensation of photochemically pro- 
duced sulfuric acid. POM and BC were directly emitted and 
accumulated to 9.0 /igm~^ and 6.0 /igm~^, respectively, un- 
til 18:00 LST when the emissions stopped. After 18:00 LST 
the mass densities declined due to dilution, especially ni- 
trate and BC for which the background concentration was 
zero. 

5.4. Aerosol distribution functions 

We take N{D) (m~*) to be the cumulative number dis- 
tribution, giving the number of particles per volume that 
have diameter less than D. Similarly, the cumulative mass 
distribution M{D) (kgm~^) gives the mass per volume of 
particles with diameter less than D, while the per-species 
cumulative mass distribution Mx{D) gives the mass per vol- 
ume of species x in particles with diameter less than D. We 
write N = N{oo), M = M(oo), and M^, = M^{oo) for the 
number, mass, and per-species mass densities, respectively. 

Given the cumulative densities, we define the number dis- 
tribution n(Z)) (m~^), mass distribution m(/3) (kgm~^) and 
per-species mass distribution mx{D) (kgm~^) by 



^ ' dlogio^ 

m{D) = — — 

dlogio-D 



m^iD) = — — ^ 
rflogio 



D' 



The initial, background, and emitted number densities used 
in this paper will all be superpositions of log-normal distri- 
butions, each defined by 



where N (m~^) is the total number density, Dgn (m) is the 
geometric mean diameter, and erg (dimensionless) is the ge- 
ometric standard deviation. 

To discuss the composition of a particle, we refer to cer- 
tain mass fractions of species, as 



/bc,pom = 

/BC.dry = 
/H2 0,all = 



AtBC + Mpom 

/JH20 

Mall ' 



(24) 
(25) 
(26) 



where we recall that Ux (m^) is the mass of species a; in a 
given particle, /iaii (m ) is the total wet mass of the particle, 
and /Udry = Mall — MH2O (m^) is the total dry mass. 

We extend the number and mass densities to be functions 
of both particle composition and diameter. That is, the 
two-dimensional cumulative number distribution Ny^x{f, D) 
(m^*) is the number of particles per volume that have a 
mass ratio of y to a; less than / and a diameter less than D. 
The two-dimensional number distribution ny,x{f,D) (m^^) 
is then defined by 



d^Ny,xif,D) 

9/aiogioD ' 



(27) 



The two-dimensional mass distribution my^x{f, D) (kgm~^) 
and cumulative mass distribution My^x{f,D) (kgm~^) are 
defined similarly. 

We can also define two-dimensional densities based on 
other particle quantities. In particular, if we denote by k 
the number of coagulation events that a given particle has 
experienced during the simulation time then we can define 
the two-dimensional singly-cumulative number distribution 
Ncoa.g{k, D) (m~^) to be the number of particles per volume 
with k coagulation events and diameter less than D. Then 
«coag(fc, -D) (m~^) is defined by 



,{k,D) 



dN,o^g(k,D) 

aiogio-D 



(28) 



For ease of comparison between different plots we fre- 
quently use normalized densities denoted by a hat, so the 
normalized two-dimensional number distribution hy^x (/, D) 
(dimensionless) is defined by 



ny,x{f, D) 
N 



(29) 
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Figure 8. Number size distributions n(D) for the sim- 
ulation with coagulation after 1, 6, and 24 hours. For 
comparison the distribution without coagulation after 24 
hours is also shown. 
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and similarly for the mass distribution. 

We also find it convenient to plot one-dimensional mass 
densities for certain composition ranges, as done in Fig- 
ure 14. We write mBC,dry{[fi, f2], D) (kgm~^) to refer to 
the total mass distribution (including water), where /ecdry 
is between /i and /2, 

mBC,dry([/l,/2],i?) = MeCdry (/2 , -D) - MeCdry (/l , I>) • 

(30) 

5.5. Aerosol size distribution evolution 

Figure 8 shows results of the number size distributions 
n{D) after 1, 6, and 24 hours of simulation including coag- 
ulation. For comparison, the result after 24 hours of simu- 
lation without coagulation is also shown. The distribution 
after 1 hour still resembled the bimodal initial distribution 
(compare Table 2), which was identical to the background 
distribution. After 6 hours the distribution was primarily 
determined by the emissions. Concurrently, condensation of 
secondary species caused aerosol growth. Particles at small 
sizes were depleted due to coagulatiou. After 24 hours the 
Aitkcn mode of the background appeared again as a result 
of dilution. Compared to the size distribution without co- 
agulation, the size distribution with coagulation showed a 
substantial decrease in number density for particles smaller 
than 0.1 fim. With coagulation the total number density N 
peaked at the end of the emission period after 12 hours with 
a maximal value of 1.67 ■ 10^" m^'\ After this. N declined 
due to coagulation and dilution to 7.36- 10'^ m^^. The simu- 
lation without coagulation lead to a maximum total number 
density of 2.39- lO"' m"'\ and a final value of 1.54- 10^" m"^. 
This means that coagulation decreased the peak and final 
total number concentrations by 30% and 52%, respectively. 
Comparing the number densities for the specific diameters 
D = 0.03, 0.05, 0.07, and 0.1 /xm with and without coagula- 
tion, we find that coagulation decreased the number density 
n{D) by 82%, 84%, 66%, and 29%, respectively 

We notice that for all size distributions shown in Figure 8 
the results are somewhat "noisy" at small and large diam- 
eters. This noise is inherent to the stochastic model that 
is used for coagulation, dilution and emissions. Towards 
the edges of the size spectrum only a few particles are being 
used to represent the size distribution due to the low number 
density. Single particle variations arising from the stochastic 
model thus appear as a noisy curve. This could be rectified 
by averaging repeated Monte Carlo simulations or by using 
a variable number of physical particles per computational 
particle, as in the Mass Flow Algorithm [Babovsky, 1999; 
Eibeck and Wagner, 2001]. 

5.6. Aerosol mixing state evolution 

While Figures 7 and 8 give an overview of aerosol size 
distribution and composition just like we obtain from trstr 
ditional distribution-based models, they do not address the 
issue of mixing state. To elucidate how the mixing state 
evolved over the course of the simulation we display the data 
as shown in Figure 10. The panels show the two-dimensional 
number distributions as a function of dry size and mass frac- 
tion of BC, /ecdry, after 1, 5, 7, and 24 hours of simulation. 
This corresponds to LST 07:00, 11:00, 13:00, and 06:00 of 
the next morning. Our definition of the two-dimensional 
size distribution is described in Section 5.4 above. 

We will discuss the evolution of the two-dimensional num- 
ber size distribution in conjunction with Figure 11. The gray 
scale in this figure shows the water content of the particles, 
/HjO.aii, as a function of BC nuxing state, /ecdry, and parti- 
cle size. We also include in Figure 9 the temporal evolution 
of the composition of three representative particles to aid 
the interpretation. These particles are labeled with PI, P2 
and P3 in Figure 10. 

Figure 10 shows the BC mixing state, /acdry, relative to 
all other dry constituents. Since even at the time of emis- 
sion no particles were pure BC, particles were not present 



at /scdry = 100%. Fresh emissions from diesel vehicles 
(/ecdry = 90%) and gasohne vehicles (/ecdry = 34%) 
appear as horizontal lines since particles in one emission 

category were all emitted with the same composition. At 
/BC.dry = 0% all the particles appear that do not contain 
any BC (i.e. background particles and particles from meat 
cooking emissions that have not undergone coagulation with 
particles contaiinng BC). After 1 hour (07:00 LST) a small 
number of particles in between these three classes indicate 
the occurrence of coagulation. 

Under the initial ambient conditions the emitted diesel 
and gasoline particles accumulated small amounts of am- 
monium sulfate, ammonium nitrate and water. After 06:44 
LST the relative humidity fell below 85%, which is the del- 
iquescence point of the inorganic mixture of ammonium, 
sulfate and nitrate. As a result of the hysteresis of parti- 
cle deliquescence and crystallization, the particles that had 
been emitted up to this point stayed wet throughout the 
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Figure 9. Time history of the composition of three in- 
dividual diesel particles, PI, P2, and P3. PI is emitted 

at 06:05 LST and always contains water. P2 is emitted 
at 06:56 LST and is initially dry but becomes wet in the 
afternoon. P3 is emitted at 16:18 LST later in the day 
when little condensation occurs. 
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Figure 10. NornicLlizcd two-dinicnsional number distribution TiBC.drv (f,D) after 1, 5, 7, and 24 hours 
of simulation. The labels PI, P2, and P3 track three individual diesel omission particles as they evolve 
over the course of the simulation. The ma:ximum plotted value for nscdry (/, is capped at 4 to allow 
better resolution. 
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Figure 11. Water fraction, /H20,aii a function of BC fraction /scdry and dry diameter after 1, 5, 7, 
and 24 hours of simulation. The labels PI, P2, and P3 track three individual diesel emission particles 

as they evolve over the course of the simulation. Note that the water fraction of wet particles is plotted 
over the hashing for dry particles and sometimes obscures it. In particular, after 1 hour (07:00 LST) 
there are dry diesel and gasoline particles present but they are not visible. The water fraction plotted 
for a given two-dimensional bin is the minimum of the water fraction for all wet particles in that bin. 
For example, after 24 hours the particle PI is very wet (see Figure 9) but there are much drier particles 
present with similar composition, giving a low /HjO.aii value on the plot at PI. The majcimum plotted 
value for /HgO.aii is capped at 50% to allow better resolution. 
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Figure 12. Normalized two-dimensional number distribution nBC,dry(/, f) after 1, 5, 7, and 24 hours 
of simulation, as in Figure 10, but for the simulation without coagulation. The maximum plotted value 
for nBC,dry(/, D) is capped at 4 to allow better resolution. 
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Figure 13. Normalized two-dimensional number distribution h^c,POM{f,D) after 24 hours of simu- 
lation (06:00 LST the following day), with and without coagulation. The maximum plotted value for 
HBC,dry(/, is Capped at 4 to allow better resolution. 



ity reached 61% in the following afternoon at 16:00 LST. 
Hence, between 06:46 and 16:00 LST, wet and dry parti- 
cles co-existed in the air parcel. Particle PI in Figures 10 
and 11 is one of the particles that was omitted early and 
stayed wet throughout the sirrmlation, whereas particle P2 
started out dry and became wet only in the afternoon. For 
the wet and dry particles different thermodynamic equilib- 
ria applied which was reflected in the different development 
of their /ec.dry values. 

As the single-particle plot Figure 9 shows, the wet par- 
ticles contained nitrate from the beginning and kept taking 
up nitrate, while during the first few hours vapor pressures 
of HNO3 and NH3 were too low to allow nitrate formation 
on dry particles. Due to this difference in nitrate formation. 



after 5 hours (11:00 LST) the wet particles appear distinct 
from the dry particles in Figure 10 and reached lower /ecdry 
values, reflecting their laiger ammonium nitrate content. 
This changed after 11:30 LST. At this time HNO3 and 

NH3 were high enough so that ammonium nitrate formed 
on the dry particles. They accumulated ammonium nitrate 
quickly, and /ecdry decreased rapidly for the dry particles. 
As a result, /ecdry of the dry particles fell below /ecdry 
of the wet particles, as is evident in the graph for 7 hours 
(13:00 LST) in Figures 10 and 11. 

After 18:00 LST the ammonium rntrate formations 
stopped, as NH3 concentration dropped to near zero (com- 
pare Figure 6). Therefore the fresh particle emissions after 
this time did not accumulate much condensable material and 
stayed at high /Bc.dry values. This is reflected in the single- 
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particle plot Figure 9, which shows the diesel particle P3 
that was emitted in the afternoon. The mass of secondary 
species for this particle was much lower than its BC content. 
After 12 hours (18:00 LST) both particle and gas emissions 
stopped, and the particle distribution changed mainly due 
to coagulation and dilution. The particle number decreased 
as a result of coagulation and continued dilution with the 
background, but this effect is not visible in the normalized 
number densities. During the evening hours the relative hu- 
midity increased again and particles took up a substantial 
amount of water. As the 24 hour plot in Figure 11 shows, 
the water content depended on the mixing state, i.e. for 
a given size we find particles with water mass-fractions be- 
tween near 0% and 66%. 
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Figure 14. Mass distribution after 24 hours (06:00 LST 
the following day) for throe different mixing state 
ranges: mBc,dry([0%, 2%], D), mBc.dry([2%, 10%], D), 
and mBC,dry([80%, 90%], D), as defined in equation (30). 
The cases with and without coagulation are plotted for 
comparison. 
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Figure 15. Two-dimensional number distribution 

ncoag(fc, D) showing the number of coagulation events ex- 
perienced after 24 hours (06:00 LST the following day), 
as defined in equation (7). The maximum plotted value 
for ncoag(fe,£') is capped at 6 • 10^ m~^ to allow better 
resolution. 



Comparing the result for the end of the simulation to 
the results at previous times, we note that at the end of 
the simulation particles below D = 0.03 /nm were depleted 
due to coagulation. A continuum of mixing states formed 
in between the extreme mixing states of /bc dry = 0% and 

/BCdry = 90%. 

The comparison to the case without coagulation gives re- 
sults as displayed in Figure 12. This figure is analogous to 
Figure 10 and shows the mixing state /scdry of BC with 
respect to the sum of all other substances. Without coagu- 
lation similar frontal features appeared, but diesel particles 
and gasoline particles remained more clearly distinct un- 
til 6 hours of simulation (12:00 LST) without the mixing 
effect of coagulation. After about 12:00 LST the mixing 
state became continuous, because the most aged diesel car 
emissions started overlapping with the relatively fresh gaso- 
line emission particles. However, since those mixing states 
were formed due to condensation only, an internal mixture 
of primary species such as POM and BC existed only when 
the particles were emitted as this mixture. After 24 hours 
of simulation without coagulation, mixed particles smaller 
than D = 0.03 /xm were still present while they were depleted 
in Figure 10 with coagulation. 

The impact of coagulation on the mixing state with re- 
spect to the primary components BC and POM is shown in 
Figure 13. The left figure displays the BC mixing state with 
respect to POM, /bc.pom, after 24 hours with coagulation. 
POM was emitted as a constituent of primary particles, 
which can be seen as horizontal lines with high number den- 
sities at /bc.pom = 90% for diesel car emissions, /bc,pom = 
34% for gasoline car emissions and /bc.pom = 0% for meat 
cooking emissions. The mixing states between these could 
only form as a result of coagulation. Since coagulation is 
most efficient between particles of different size, we observe 
that these mixed particles preferentially formed in a specific 
size range. For sizes larger than D = 0.05 /xm, POM/BC 
mixtures of various degree of mixing formed due to coagula- 
tion. Below 0.05 fim, coagulation produced very few parti- 
cles, so particles at these sizes were at their initial BC/POM 
mixing state. Nearly all particles below 0.03 /xm were re- 
moved by coagulation. For comparison the right panel of 
Figure 13 shows the BC mixing state with respect to POM, 
/bc.pom, at the end of the simulation without coagulation. 
For this case the intermediate mixing states did not occur, 
and particles below D = 0.03 fim remained. 

Figure 14 shows the one-dimensional size distributions of 
total mass density (including water) for different ranges of 
mixing states at the end of the simulation, comparing the 
cases with and without coagulation. From this we sec that 
coagulation did not simply reduce the number densities, but 
also shifted black carbon mass within the diameter- /bc, dry 
space. The mass of particles smaller than D = 0.05 /im 
with high BC content (/BC.dry between 80% and 90%) was 
reduced by 90% due to coagulation. The mass of particles 
smaller than D = 0.05 /iin with very low BC content (/BC.dry 
between 0% and 2% BC) was reduced by 81% when coagula- 
tion was included. A very large difference between the cases 
with and without coagulation occurred for /BC.dry between 
2% and 10% and for the size range above 73 = 0.1 /im. Mass 
in this range of parameters arose mainly from coagulation 
of large, BC-free particles with small BC-containing parti- 
cles and this mass increased by 1458% when coagulation was 
included. 

With PartMC it is straightforward to track the number of 
coagulation events experienced by the individual particles. 
Figure 15 shows the two-dimensional number distribution 
n(fc, D), with k being the number of coagulation events. At 
the end of the simulation, 5% of particles had undergone at 
least 5 coagulation events. The largest rmmber of coagula- 
tion events was 19 and occurred for particles on the right 
hand side of the spectrum with sizes larger than D = 0.4 /xm. 
Given a certain size, the number of coagulation events var- 
ied, which shows the stochastic nature of the coagulation 
process. The range of variation was greater for larger par- 
ticles. For example, while the number of coagulation events 
varied between and 5 for a D = 0.1 particle, it ranged 
between and 12 for a D = 0.3 fjxn particle. 
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6. SumniEiry 

In this paper wc presented the development and the 
apphcatiori of a stochastic particle-resolved aerosol model, 
PartMC-MOSAIC. It explicitly resolves the composition 
of individual particles in a given population of different 
types of aerosol particles, and accurately tracks their evo- 
lution due to emission, dilution, condensation and co- 
agulation. To make the direct stochastic particle-based 
method practical, we implemented an accelerated coagu- 
lation method. With this method we achieved optimal 
efHciency for applications when the coagulation kernel is 
highly non-uniform, as is the case for many realistic en- 
vironments. The highly accurate treatment of aerosol dy- 
namics and chemistry makes PartMC-MOSAIC suitable for 
use as a numerical benchmark of mixing state for more 
approximate models. The current version of PartMC is 
available under the GNU General Public License (GPL) 
at littp : / / www. mechsc .uiuc.edu/ research / mwest / part mc / , 
and the MOSAIC code is available upon request from 
R. A. Zaveri. 

PartMC-MOSAIC was applied to an idealized example 
urban plume case to simulate the evolution of urban aerosols 
of different types due to coagulation and condensation, fo- 
cusing on the aging process of BC. For the first time results 
of the aerosol composition and size distribution are avail- 
able as a fully multi-dimensional size distribution without 
any a priori assumptions about the mixing state. This detail 
of information was only achievable with a particle-resolved 
model. 

To display the results, we projected the multi-dimensional 
size distributions to two-dimensional distributions depend- 
ing on particle size and BC mass ratio. We specifically 
discussed the results for BC mass ratios defined with re- 
spect to all other dry constituents, /scdry, and to POM, 
/bc.pom- Due to the diurnal variations in temperature, rel- 
ative humidity, and gas phase concentrations, the thermo- 
dynamic equilibrium conditions for the ammonium-sulfate- 
nitrate system changed continuously. The aerosol hydration 
hysteresis effect led to the co-existence of metastable (wet) 
and stable (dry) particles in the air parcel during the day- 
time, depending on their time of emission. Since the for- 
mation of ammonium nitrate depends on the particle phase 
state, this in turn resulted in pronounced differences in how 
the aging proceeded. As a result of coagulation and conden- 
sation, after 24 hours of simulation, the aerosol population 
evolved into a state where a continuum of BC mixing states 
existed. Coagulation was effective in removing smaller par- 
ticles, reducing number densities by 84%, 66%, and 29% at 
diameters of 0.05, 0.07, and 0.10 //m in the example simula- 
tion. 
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